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We study the finite temperature transition in QCD with three flavors of equal masses using the R 
^"*^ \ and RHMC algorithm on lattices with temporal extent N T = 4 and 6. For the transition temperature 

in the continuum limit we find roT c — 0.429(8) for the light pseudo-scalar mass corresponding to 
the end point of the 1st order transition region. When comparing the results obtained with the 
R and RHMC algorithms for p4fat3 action we see no significant step-size errors down to a lightest 
pseudo-scalar mass of m pa ro — 0.4. 
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I. INTRODUCTION 



Lattice QCD has established the existence of a transition from hadron gas to a new state of strongly interacting 
matter where quarks and gluons are no longer confined inside hadrons and which is usually called the quark gluon 
plasma 0, The nature of this transition depends on the quark content and quark masses. For infinite or very 
large quark masses the transition is a 1st order deconfining transition. In the opposite case of zero quark masses 
one may have a 2nd order chiral phase transition for 2 flavors or a 1st order chiral phase transition for 3 flavors. 
, For intermediate masses the transition is just a rapid crossover, meaning that thermodynamic quantities change very 
rapidly in a narrow temperature interval. The boundary of the 1st order transition region of 3 flavor QCD as a 
function of mass has been studied using improved (p4) and standard staggered actions [J, 0] . There is a significant 
■ discrepancy regarding the value of the quark masses, or equivalently the value of the pseudo-scalar meson masses where 
the transition changes from crossover to 1st order. With an improved action it was found that the 1st order transition 
ends for a pseudo-scalar meson mass of about 70 MeV, while with the standard action it ends for pseudo-scalar meson 
masses of about 190 MeV or larger Q. 

It has been observed that the pressure and energy density normalized by its ideal gas value shows almost the 
same behavior as a function of T/T c for SU(3) pure gauge theory, 2 flavor, 2+1 flavor and 3 flavor QCD [l|. Thus 
flavor and quark mass dependence of these quantities in the first approximation, is determined by flavor and quark 
mass dependence of transition temperature T c . Therefore it is very instructive to study the flavor dependence of the 
transition temperature. Such a study has been performed in Ref. 0] on lattices with temporal extent N T = 4 using 
the improved staggered p4 action. 

In the past most simulations with 2 and 3 flavors of staggered fermions have been done using the hybrid molecular 
dynamics R (HMDR) algorithm [j§], often called simply the R-algorithm. It has finite step-size errors of (D(dt 2 ) where a 
step-size dt is used in the molecular dynamics evolution. Recently the rational hybrid Monte Carlo (RHMC) algorithm 
has been invented which allows simulations of theories with fractional powers of the fermion determinant, for example 
2 and 3 flavors of staggered fermions, without finite step-size errors [9j. Therefore the most recent thermodynamics 
studies use the RHMC algorithm [l(| [ll|, [H, 03 ■ R has been observed in Ref. [ll[ that the use of the exact RHMC 
algorithm reduces the value of the critical quark mass where the transition turns to 1st order by 25% in the case of 
the standard staggered action. 

The purpose of this paper is twofold. First we would like to study the transition in 3 flavor QCD, extending 
the previous studies to smaller quark masses and smaller lattice spacings using the improved p4 staggered fermion 
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action. Second we would like to compare the R-algorithm with the new RHMC algorithm. The rest of the paper is 
organized as follows. In section II we discuss the calculational set-up. In section III we show our results for finite 
temperature calculations. Section IV discusses the zero temperature simulations needed to set the scale and the 
transition temperature in physical units. In section V we present a comparison of the R and RHMC algorithms. 
Finally, section VI contains our conclusions. A technical discussion of different fat link actions is given in Appendix 
A. In Appendix B we discuss the eigenvalues and eigenvectors of the p4 action in the free field limit. 

II. LATTICE FORMULATION AND SETUP 

Most of the simulations discussed in this paper were done using the p4fat3 action, which is also simply called the 
p4 action • To improve the rotational symmetry, which is violated on the lattice, bent 3- link terms are added to the 
1-link term of the standard staggered action [l5j]. Properly chosen coefficients of the 1-link and the 3-link terms can 
eliminate, at tree level, the 0(a 2 ) errors in the dispersion relation for staggered fermions The violation of flavor 
symmetry which is present in the staggered fermion formulation can be significantly reduced by replacing the normal 
link in the 1-link term by a fat link which is a sum a of normal link and 3-link staples [l6l |. This 3-link fattening is 
the origin of the name p4fat3. Though this type of fat link action is a big improvement over the standard staggered 
fermion action, further improvement of the flavor symmetry can be obtained by adding 5- and 7-link staples [l^ |. 
In fact one can eliminate the effect of flavor symmetry breaking at order 0(g 2 a 2 ) using a suitable combination of 
3-, 5- and 7- link staples leading to what is called the fat7 action [17J. We have also done calculations with the p4 
action with fat7 fat links. Unfortunately it turns out that on the coarse lattices used in our study of 3 flavor QCD 
thermodynamics this action has some undesirable features. It leads to the occurrence of a bulk transition, which we 
will discuss in more detail in Appendix A. 

To study staggered fermions with less than four flavors, we use the rooting procedure, i.e. each fermion flavor is 
represented by (detAd) 1 ^ 4 , where M. is the staggered fermion Dirac operator. For a recent discussion of this procedure 
see Ref. [l8| . Most of our simulations have been done using the standard R-algorithm Q . As in Ref. Q the step-size 
of the molecular dynamics evolution was set to dt = m/2.5 for the staggered quark mass m. 

We also performed calculations with the RHMC algorithm. In this algorithm an optimal rational approximation 
is used to evaluate the fractional power of the determinant [9(. More precisely one finds the optimal approximation 
for (M.^ M) v where M. = 2m + D is the usual staggered fermion matrix and for the three flavors v — 3/8. Using 
sufficiently high order polynomials the rational approximation can be made arbitrarily precise for the given spectral 
range of the fermion operator. For the standard staggered fermion action the spectral range of M. is well known, 
the smallest eigenvalue of this matrix is 4m 2 . The largest eigenvalue can be estimated in the free field limit to be 
^max = 16 + 4m 2 . For the case of the p4 staggered action the smallest eigenvalue is the same, while the largest 
eigenvalue in the free case is A 2 na;r = 50/9 + 4m 2 . In appendix B we give the derivation of this result. It turns out that 
in all our simulations the largest eigenvalue was smaller than 5.0, so we choose A 2 na2 . = 5.0 as the upper limit on it. 
For the range of the quark masses studied by us, which include quark masses as light as l/20th of the strange quark 
mass, it is sufficient to use polynomials of degree 12 to achieve machine precision with the rational approximation. 
Therefore in the Metropolis accept /reject step we used polynomials of degree 12. In the molecular dynamics evolution 
we used a less stringent approximation of the determinant since any errors in the evolution, including the dt 2 step-size 
errors, are eliminated by the accept/reject step. It has been found that in most cases it is sufficient to use polynomials 
of degree 5-6 without compromising the acceptance rate. Furthermore, the stopping criteria for the conjugate gradient 
inversions can be relaxed to 10~ 5 in the molecular dynamics evolution without significant effect on the acceptance 
rate. In the Monte-Carlo accept/reject step we typically use 10~ 12 for the conjugate gradient stopping criteria. The 
length of the trajectory was tmd — 0.5 in units of molecular dynamics time. 

The gauge fields and quarks contribute different amounts to the force in the molecular dynamics evolution. The 
contribution of the gauge fields is larger than that of quarks. On the other hand the cost of the evaluation of the 
force coming from the gauge fields is small. The opposite is true for the part of the force coming from the fermions. 
Therefore it is reasonable to integrate the gauge force on finer molecular dynamics time scale than the fermion force 
(l9| . In our simulations we typically used 10 gauge field updates per fermion update. The 3 flavors of staggered 
fermions are simulated as 1+1+1, i.e. we used a factor (detM.* 1 A4) 1 / 8 for each fermion flavor. Although this increases 
the number of inversions, the reduced force allows for a larger step-size dt for the same acceptance rate. The step-size 
dt was chosen such that the acceptance is about 70%. To achieve this, the stepsize was typically of the order of the 
strange quark mass used in our 2+1 flavor study [131 ]. 
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m 


N„ 


# j3 values 


max. no. of traj. 


4 


0.100 


16 


4 


42000 




0.050 


8, 16 


8, 14 


4130, 2650 




0.025 


8, 16 


8, 9 


6250, 8650 




0.010 


8, 16 


9, 6 


2460, 4660 




0.005 


12, 16 


11, 8 


1360, 3000 


6 


0.100 


16 


16 


6000 




0.050 


16 


14 


7900 




0.020 


16 


10 


16000 




0.010 


16 


9 


10900 



TABLE I: Parameters of the numerical simulations 



III. FINITE TEMPERATURE SIMULATIONS 



Our studies of the QCD transition at finite temperature have been performed on lattices of size N% x N T . The 
lattice spacing, a, relates the spatial (N a ) and temporal (N T ) size of the lattice to the physical volume V = (N a a) 3 
and temperature T — l/N T a, respectively. The lattice spacing, and thus the temperature, is controlled by the gauge 
coupling, j3 — 6/<? 2 , as well as the bare quark masses. The parameters of our finite temperature simulations are given 
in Table U We extended the results of Ref. Q in two respects. Compared to Ref. (tJ we have added a smaller mass 
value to = 0.005 for N T — 4 and extended the runs at larger quark masses to achieve a better statistical accuracy. In 
addition we have studied the finite temperature transition on N T = 6 lattices for two values of the quark mass to. All 
the results presented in this section have been obtained with the R algorithm. 

As mentioned in the introduction, in 3 flavor QCD for small quark mass we have a 1st order phase transition which 
turns into rapid crossover at the quark mass corresponding to the light pseudo-scalar mass of about 70 MeV [J] . The 
transition is signaled by a rapid change in bulk thermodynamic observables (energy density pressure) as well as in 
the chiral condensates and the Polyakov loop expectation value, 

w = (^ Tr £ ft two) ■ ( 2 ) 

\° iV ^ x x = l I 

which are order parameters for a true phase transition in the zero and infinite quark mass limit, respectively. Note 
that we have defined the chiral condensate per flavor degree of freedom, hence the factor 1/3 in Eq. JTJ) . 

We use the Polyakov loop susceptibility as well as the disconnected part of the chiral susceptibility to locate the 
transition temperature. 

XL = iV 3 ((L 2 }-(L) 2 ) , (3) 

We calculate the value of the Polyakov loop at the end of every trajectory. For each tenth trajectory we calculate 
the value of 'tpip an d Xq using ten Gaussian random vectors. In Fig. [1] we show the disconnected part of the chiral 
susceptibility calculated on our N T — 4 lattice with different spatial volumes at different quark masses. In Fig. [5] we 
show the Polyakov loop \l and chiral \q susceptibilities calculated on 16 3 x 6 lattice. The location of peaks in the 
susceptibilities has been determined using Ferrenberg-Swedsen re-weighting for several values (3 in the vicinity of the 
transition. Errors on the peak location have been obtained from a jackknife analysis where Ferrenberg-Swedsen re- 
weighting has been performed on different sub-samples. The resulting pseudo-critical couplings are shown in Table HT1 
In finite volume the pseudo-critical couplings (3 C determined from the Polyakov loop correlator and chiral susceptibility 
are generally different. In the case of the crossover this difference can persist even in the infinite volume limit. From 
Table HIl we see that in most cases the two pseudo-critical couplings are identical within statistical errors even for small 
volume 8 3 x 4. The cases where this difference is the largest are the cases where (3 c<q has large statistical errors. For 
example for the 16 3 x 6 lattice and to = 0.05 we find (3 c .l — (3 c .q — .0113(328). Therefore we have also calculated the 
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3.25 3.30 3.35 3.40 3.25 3.30 3.35 3.40 

FIG. 1: The disconnected part of the chiral susceptibility calculated on N T = 4 lattices at different quark masses. 



60 




3.40 3.45 3.50 3.55 3.60 3.65 3.70 3.40 3.45 3.50 3.55 3.60 3.65 3.70 

FIG. 2: The disconnected part of the chiral susceptibility (left) and the Polyakov loop susceptibility (right) calculated on 16 3 x 6 
lattices. 

weighted average of j3 c ,L and (3 c ,q which is shown in the last column of Table ITT1 together with the corresponding error. 
This error is calculated from the statistical errors and the difference between the central values added quadratically. 
The difference in pseudo-critical couplings determined on 8 3 x 4 and 16 3 x 4 lattices is typically small, indicating 
small finite volume effects. As a general tendency the pseudo-critical coupling (3 C shifts toward smaller values with 
increasing volume. In agreement with earlier calculations we find that the position of peaks in Xq arL d XL show only 
little volume dependence and that the peak height changes only little, although the maxima become somewhat more 
pronounced on the larger lattices. This is consistent with the transition being a crossover rather than a true phase 
transition in the infinite volume limit for the range of quark masses explored by us. 
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m 


iV CT 


(3c,l [from xl] 


Pc, g [from Xq] 


/3c [averaged] 


4 


0.100 


16 


3.4800(27) 


3.4804(24) 


3.4802(18) 




0.050 


16 


3.3884(32) 


3.3862(47) 


3.3877(34) 






8 


3.4018(35) 


3.3930(201) 


3.4015(94) 




0.025 


8 


3.3294(27) 


3.3270(28) 


3.3283(31) 




0.010 


16 


3.2781(7) 


3.2781(4) 


3.2781(3) 






8 


3.2858(71) 


3.2820(61) 


3.2836(60) 




0.005 


16 


3.2656(13) 


3.2678(12) 


3.2667(24) 






12 


3.2659(13) 


3.2653(12) 


3.2656(10) 


6 


0.200 


16 


3.8495(11) 


3.9015(279) 


3.8495(520) 




0.100 


16 


3.6632(55) 


3.6855(105) 


3.6680(228) 




0.050 


16 


3.6076(24) 


3.6189(328) 


3.6077(115) 




0.020 


16 


3.4800(110) 


3.4800(80) 


3.4800(65) 




0.010 


16 


3.4518(50) 


3.4510(83) 


3.4516(44) 



TABLE II: Critical couplings determined from the location of peaks in the Polyakov loop susceptibility as well as in the 
disconnected parts of the chiral susceptibilities. The last column gives the average of j3 c ,L and /3 c ,q with combined statistical 
and systematic errors. 



IV. ZERO TEMPERATURE CALCULATIONS AND THE TRANSITION TEMPERATURE 

In order to determine the transition temperature in units of some physical quantity we performed zero temperature 
calculations on 16 3 x 32 lattices in the vicinity of the pseudo-critical coupling (3 C . The parameters of these calculations 
together with the accumulated statistics are summarized in Table IIIII We have calculated the static quark potential 
and meson correlators on each 10th trajectory generated. 

The static potential has been calculated using the ratios of the Wilson loops at two neighboring time-slices and 
extrapolating them to infinite time separation with the help of constant plus exponential form. The spatial transporters 
in the Wilson loop have been constructed from spatially smeared links with APE smearing. The weight of the 3 link 
staple was 7 = 0.4 and we used ten steps of APE smearing. From the static potential we have determined the string 
tension and the Sommer parameter ro defined as [20] 



dr 



= 1.65. (5) 



-To 



When extracting ro and the string tension on coarse lattices, such as the ones used in thermodynamics studies, the 
violation of rotational symmetry has to be taken into account. We do this using the procedure described in detail 
in our recent paper jl3[ . The value of Sommer scale and the string tension are given in Table IIIII for different quark 
masses. Having determined ro for different gauge couplings and quark masses allows us to perform interpolations of 
rp/a in these parameters. As in Ref. [L3( ] we use the following renormalization group inspired interpolation ansatz 

(ro/a)- 1 = R(/3){1 + Ba 2 {(3) + Ca i {(3))e A( - 2m ' +m ^ +D . (6) 

Here R(f3) is 2-loop beta function of 3 flavor QCD. In the interpolation we also used the values of ro/a determined in 
our 2+1 flavor study [H in addition to those shown in Table UTIl giving A = 1.45(5), B = 1.20(17), C = 0.21(6) and 
D = 2.41(5) with x 2 /dof = 0.9. This was the reason for using the notation mi and m s for the light and the strange 
quark mass in Eq. ([5]) . For 3 degenerate flavors of course m s = mi = m. In Table IIIII we also show the values of 
ro/a obtained from this ansatz for each of the parameter sets. 

Meson masses have been calculated using the four local staggered meson operators. We used point-wall meson 
correlation functions with a Z2 wall source. To extract the meson masses from the correlation functions we used 
a double exponential ansatz which takes into account the two lowest states with opposite parity. The two lowest 
pseudo-scalar meson masses as well as the lightest vector meson mass are shown in Table IIIII The breaking of the 
flavor symmetry can be quantified by the quadratic splitting of the pseudo-scalar masses: A ps — {m 2 s2 — 'm<p S )'r 2 , 
where m pS 2 is the mass of the lightest pseudo-scalar non-Goldstone meson that is present with staggered fermions. 
This quantity should be quark mass independent for sufficiently small quark masses and should vanish as 0(a 2 ) when 
the continuum limit ( a — > 0) is approached. In the last column of Table [Till we show the value of A ps from our scale 
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m 


# traj 


trips 


m ps 2 


my 


ro 


(^*o) smooth 






3.3877 
3.3270 
3.2680 


0.050 
0.025 
0.005 


7800 
12000 
1500 


0.7084(1) 
0.5118(3) 
0.2341(9) 


1.094(7) 
0.998(24) 
0.860(90) 


1.310(20) 
1.222(32) 
1.250(50) 


2.066(7) [7] 
1.982(14)[13] 
1.888(15) [9] 


2.061 
1.989 
1.888 


0.552(12)[12] 
0.564(11)[11] 
0.587(17)[17] 


2.97(7) 
2.90(20) 
2.44(55) 


3.46345 
3.4400 


0.020 
0.010 


4420 
4290 


0.4413(8) 
0.3210(7) 


0.665(5) 
0.594(7) 


0.908(11) 
0.882(20) 


2.797(20) [20] 
2.770(13) [13] 


2.813 
2.779 


0.404(6) [6] 
0.405(6) [6] 


1.94(9) 
1.92(7) 



TABLE III: Parameters of the zero temperature simulations, meson masses, the Sommer scale ro and the string tension. Also 
shown is the value of ro obtained from the interpolation formula ©. The upper part of the table refers to scale setting runs 
for our N T = 4 lattices while the lower part to our N T = 6 calculations. In the last column the splitting between the lightest 
non-Goldstone and the Goldstone pseudo-scalar meson masses squared is shown. All dimensionfull quantities are given in units 
of the lattice spacings. 



setting run for N T — 4 and N T = 6. As we see from the table this quantity does not decrease quite as fast as a 2 . This 
is an indication that on the coarse lattice, corrections to asymptotic scaling arc still important. 

With the help of the interpolation formula we can calculate the lattice spacing in units of ro and thus r T c for 
different pseudo-scalar meson masses, which is shown in Fig. The error in a(f3 c ) results in an error in the value 
of tqT c which is shown in Fig. [3] as thin error-bars. The uncertainty in (3 C itself also contribute to the uncertainty in 
tqT c1 which is shown as a thick error-bar in the figure. 

If there is a critical point in the (to, T)-plane, then universality dictates that T c (m) — T c (m e ) ~ (m — rrf) 1 ^ 5 ^ 
with fj and 6 being critical exponents. In the case of three degenerate flavors the line of the 1st order transition in the 
(to, T) plane ends in a critical end-point m e ,T c (m e ) belonging to the Z(2) universality class. For this universality class 
we have 8(3 = 1.5654. Therefore we attempted a combined continuum and chiral extrapolation using the following 
extrapolation ansatz 

r T c (m ps ,N T ) = r T c \ cont (m ps ) + A ((r m ps ) - (r TO pS:C ) J + B/N r . (7) 

The value of the quark mass where the transition changes from 1st order to crossover, i.e. the mass corresponding to 
the end-point, has been estimated in 0] using N T — 4 lattices to be m e — 0.0007(4). This translates into the value of 
the pseudo-scalar mass 

r TO; s =0.16l 3 5 . (8) 

It turns out that this large uncertainty in the value of m ps produces an uncertainty in the extrapolated value of 
T c (m ps ) which is much smaller than the statistical errors. The extrapolation according to Eq. iJJJ) yields 

T (m e ) 

r T c (m; s ) = 0.429(8), cV =0.391(9). (9) 

For the fit to tqT c data we get x 2 /dof = 0.7, while for the fit to T c {mtg) / ' \Jo data we have x 2 /dof — 0.4. The quark 
mass dependence of the transition temperature is described by Eq. |J7J) only for m > m e . For smaller quark masses 
the transition is first order and T c depends linearly on the quark mass, i.e. we expect T c (m ps , N T ) — T c \^^ t al + 
Arrip S + B/N 2 . If we would insist on the linear dependence of the transition temperature on the quark mass in the 
entire mass range, the combined chiral and continuum extrapolation would give 

r T c = 0.419(9), = 0.383(10). (10) 

The x 2 /dof we have found for this fit is almost the same as for the one above. The value of T c /y/a is slightly 
smaller than the estimate of Ref. [7| based on N T — 4 lattice and larger quark masses. This is due to the continuum 
extrapolation performed in the present work. Note, however, that the value T c (0) * =4 / \fa = 0.417(9) is entirely 
consistent with Ref. Q- The value of T c could be compared with the corresponding 2 + 1 flavor value T c tq = 
0.444(6) [+12] [—6] in the limit of vanishing u and d quark masses but fixed physical value of m s [l3|. Thus the flavor 
dependence of r$T c is about or smaller than 5%. One should also note that the difference between the transition 
temperature calculated on iV r = 4 and N T = 6 lattices is very similar to that found in 2+1 flavor case [l3|. In Ref. 
the transition temperature in the chiral limit has also been estimated in units of the vector mass. Our estimate 
for T c /mv\m=o is consistent with that result. 
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0.50 
0.48 
0.46 
0.44 
0.42 
0.40 
0.38 



0.36 



T c /Vo 




0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 



0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 



FIG. 3: The value of the transition temperature in units of ro (left) and in units of <Jd (right) calculated on iV r = 4 and N T = 6 
lattices as function of m ps together with continuum extrapolated value. The vertical line and band indicates the value of m ps 
where the transition becomes 1st order. 




FIG. 4: The comparison of the disconnected part of chiral condensate susceptibility calculated with the R and RHMC algorithms 
for p4fat3 with m = 0.01. 



V. COMPARISON OF R AND RHMC ALGORITHM 



We investigated the effect of the finite step-size errors of the R algorithm on the properties of the finite temperature 
transition. We performed calculations on 8 3 x 4 lattices using the p4fat3 action as well as the p4fat7 action and the 
later will be described in Appendix A in more detail. In the calculations with the p4fat3 action we used a quark mass 
of to = 0.01, while in case of p4fat7 action we used two quark masses m = 0.1 and 0.035. In our calculations with the 
R algorithm the step-size of the molecular dynamics evolution was set to be dt = m/2.5. Some additional calculations 
have been done at twice smaller step-size dt = to/5. We have calculated the chiral condensate and the Polyakov loop 
and determined the pseudo-critical couplings which are summarized in Table IIVI In Fig. [4] we compare the chiral 
condensate susceptibility calculated using the R-algorithm and RHMC algorithm for the p4fat3 action. We find that 
for the p4fat3 action the results obtained with R and RHMC algorithms are identical within statistical errors. 

The situation is different for p4fat7. In Fig. [5]the expectation value of the Polyakov loop and the chiral condensate 
calculated with the two algorithms are shown for two values of the quark mass. Here we see significant differences in 
the value of the chiral condensate and Polyakov loops calculated with the R algorithm and step-size dt = m/2.5 and 
the corresponding result obtained with RHMC algorithm. We see also a small but statistically significant difference 
in the value of the pseudo-critical coupling calculated with the two algorithms, c.f. Table ITVl The difference becomes 
much less visible when the step-size is decreased to dt = to/5. Figure [5] also suggests that the difference between the 
results of the two algorithms becomes larger for the smaller quark mass. 

For the p4fat7 action we also performed a zero temperature calculation on 16 3 x 32 lattices using the RHMC 
algorithm to determine the scale. We have calculated meson masses as well the static quark potential. From the 
later we have extracted ro . The results of these calculations are also given in Table IIVI Now we can estimate the 
transition temperature in units of ro for to = 0.035 calculated with the two algorithms. For the R-algorithm we get 
r T c = 0.542(3) while for the RHMC algorithm we have r T c = 0.552(3). Thus in the case of the p4fat7 action the R 
algorithm underestimates the transition temperature roughly by 2%. 



8 




FIG. 5: The Polyakov loops (upper part) and the chiral condensate (lower part) calculated with RHMC and R algorithms for 
the p4fat7 action with m = 0.1 and m = 0.035. 



Action 


m 


Algorithm 






trips 


ro /a 


p4fat3 


0.010 


HMDR 
RHMC 


3.2858(71) 
3.2820(11) 


3.2820(61) 
3.2820(11) 






p4fat7 


0.100 


HMDR 
RHMC 


2.9850(25) 
2.9939(33) 


2.9753(53) 
2.9831(20) 






p4fat7 


0.035 


HMDR 
RHMC 


2.7514(6) 
2.7540(6) 


2.7485(7) 
2.7515(7) 


0.7884(5) 
0.7897(7) 


2.1661(123) 
2.2063(108) 



TABLE IV: Comparison of the R and the RHMC algorithms for the pseudo-critical couplings and the scale at the transition 
temperature. 



VI. CONCLUSIONS 



In this paper we have studied the phase transition in 3 flavor QCD at finite temperature using N T = 4 and N T = 6 
lattices. For the quark mass corresponding to the second order end-point we find the critical temperature to be 
roT c = 0.429(8). The transition temperature in the chiral limit is about 2% smaller than the above value. For a given 
pseudo-scalar meson mass the difference between the transition temperature in 3 flavor and 2+1 flavor case is less 
than 5%. We also find that the cut-off dependence of the transition temperature in 3 and 2+1 flavor QCD is very 
similar. Furthermore, we find that finite step-size errors present in the R algorithm are negligible, at least for the 
p4fat3 action at the quark masses studied. 
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FIG. 6: The chiral condensate calculated on N T = 4 (left) and N T = 6 (right) lattices for p4fat7 action. 
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Appendix A 

In this appendix we are going to discuss the properties of the finite temperature transition in the case of the p4fat7 
action and compare it to the p4fat3 case. In general the gauge transporter in the 1-link term of the p4 action can 
be replaced by combination of the link variable and different staples, called the fat link, without changing the naive 
continuum limit. This is true provided the coefficient of different terms in the fat link satisfy appropriate normalization 
conditions. For example in the case of the fat link with the three link staple only, this condition reads c\ + 6C3 = 1. 
Introducing five and seven link staples in addition to the three link staples give the so-called fat7 link [ljj • In this case 
the normalization condition reads: c\ + 6C3 + 24cs + 48c7 = 1. It is possible to eliminate the leading order coupling to 
the high momentum gluons with momenta (0, 7r, 0, 0), (0, 7r, 7r, 0) and (0, 7r, 7r, 7r), i.e. to suppress the flavor changing 
interaction at order g 2 a 2 if the coefficients are chosen as [ITT ] 

C-i 1 Ck 1 C7 1 . 

— = -, — = -, — = — ■ 11 

ci 2' ci 8' ci 48 K ' 

This gives then the value c\ = — 1/8 for the coefficient of the 1- link term to get the naive continuum limit. 

In Fig. [6] we show the chiral condensate calculated on the N T = 4 and N T = 6 lattices. We can see that for small 
quark masses the transition becomes strongly first order. The value of the chiral condensate in the low temperature 
phase is also much larger than for the p4fat3 action discussed in the main text. For the smallest quark mass the 
discontinuity in the chiral condensate is about the same for N T = 4 and N T — 6 lattice, although we would expect 
it to decrease by roughly a factor of (6/4) 3 ~ 3 when going from N T = 4 to N T — 6. This could mean that we are 
dealing with a bulk transition. In Fig. [7] we also compare the pseudo-critical couplings for the p4fat3 and p4fat7 
actions. We see that for large quark masses the N T dependence of the pseudo-critical coupling is similar, though their 
values are significantly different. For small quark masses the pseudo-critical couplings calculated for N T = 4 and 6 
come very close together, again suggesting that the transition may be a bulk transition. We did calculations also on 
8 4 lattice at m = 0.01. In Fig. [7]we compare the chiral condensate calculated on 8 3 x 4, 16 3 x 6 and 8 4 lattices. We 
see a sharp drop in the value of the chiral condensate, which occurs at the same /3 C for N T — 6 and 8. This again 
indicates a bulk transition. 

One may wonder which feature of the p4fat7 action is responsible for the bulk transition. The main difference of the 
p4fat7 action compared to p4fat3 action as well as to other fat link action (e.g. ASQTAD) is the negative sign of the 
one link term. Close to the continuum limit the normalization condition c\ + 6C3 + 24cs + 48c7 = 1 should insure that 
the combination of 1-, 3-, 5- and 7-link terms will describe a conventionally normalized, positive Dirac kinetic energy. 
However, at stronger coupling where the gauge field are more disordered, the staples with many links are expected 
to give a significantly smaller contribution and the 1-link term may dominate, resulting in an effective kinetic energy 
term with a possibly negative sign. While this would simply correspond to non-standard sign and normalization 
conventions for the Dirac kinetic energy, it raises the possibility that this effective kinetic energy term will change sign 
as one passes from strong to weak coupling. Such a sign change could induce a bulk transition. In addition, the change 
in magnitude of the coefficient in the effective kinetic energy (small for strong coupling and large for weak coupling) 
would appear reversed in the chiral condensate (large for strong coupling and small for weak coupling) consistent with 
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FIG. 7: The pseudo-critical couplings for the p4fat3 and p4fat7 actions (left) and the chiral condensate calculated with the 
p4fat7 action on 8 3 x 4, 16 3 x 6 and 8 4 lattices at m = 0.01 (right). 




FIG. 8: The chiral condensate (left) and its susceptibility (right) for the p4fat7' action calculated on 8 3 x 4 lattice at m = 0.01. 
We also show the chiral condensate calculated with p4fat7 and p4fat3 actions at the same quark mass. The data points for 
p4fat3 action have been shifted horizontally by A/3 — —0.437 for better visibility. The band in the right figure corresponds to 
Ferrenberg Swendsen re- weighting. 



the observed behavior. To verify this we did calculations with p4fat7 action but with different coefficients which we 
call the p4fat7' action. The coefficients were chosen to be 

— £ i £i — i — 1 £l — J_ (i9\ 
Cl ~ 4 ' 8' a. ~ 2' ci " 8' a ~ 48' 1 ' 

For this action we found no evidence for a strong first order transition but only a crossover. This can be seen for 
example in the behavior of the chiral condensate shown in Fig. [8] Both the value of the chiral condensate and the 
location of the transition point is very similar to that of the p4fat3 action. 

We also calculated the eigenvalues A = iX' + 2m of the p4fat7 Dirac operator. The normalized distribution of the 
lowest 50 eigenvalues A' is shown in Figs. l9l and [T0l for the quark masses m = 0.1 and m = 0.01, respectively. We 
have used 100 configurations for m — 0.01, and 200 configurations for m = 0.1. Given the above definition of A' the 
breaking of the chiral symmetry manifests itself in a non-zero density at A' ~ 0. In Fig. [9] we show the eigenvalue 
distribution for the larger quark mass m — 0.1. It shows the expected features: large density of eigenvalues near 
A ~ in the low temperature phase {(3 — 2.96), significant drop of eigenvalue density around zero at the transition 
(f3 = 3.0) and zero density of eigenvalues at the origin for the deconfincd phase ((3 — 3.04). The situation is different 
for the smallest quark mass m = 0.01, where we see non-zero density of eigenvalues at A ~ even in the deconfincd 
phase (P = 2.66). The large decrease in the density of eigenvalues at the origin when going from the confined phase 
(/3 = 2.635) to the deconfined explains the large drop in the value of the chiral condensate. 
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FIG. 9: The eigenvalue distributions of the p4fat7 Dirac operator for m — 0.1 calculated on 8' x 4 lattice in the confined phase 
(/3 = 2.96), at the transition (/3 = 3.00) and in the deconfined phase (/3 = 3.04). 
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FIG. 10: The eigenvalue distributions of the p4fat7 Dirac operator for m = 0.01 calculated on 8 X 4 lattice in the confined 
phase (/3 = 2.635), at the transition (/3 = 2.64) and in the deconfined phase (/3 = 2.66). 



Appendix B 

In this appendix we discuss the calculations of the largest and the smallest eigenvalue of the staggered Dirac 
operator in the free field limit. Let us start our discussion with case of the standard staggered fermions. The free- field 
staggered Dirac operator acting on a single-component fcrmion field is given by 



Mx(x) = ^2 ti m (x) (x(x + n) - x(x - M)) + 2m x( x ), 



(13) 



where rj^x) = (— l) 3! o+ x i+— x f-i are the normal staggered phases. Consider M. acting on a momentum eigenstate 

X(x) = X (p)e lp - X , (14) 
Mx(x) = J2 [x(p) e ip ' x (e lp ^ -e^Je^- 1 ] + 2m X (p) e lp ' x , (15) 



where ir = (0, 0, 0, 0), it\ = (n, 0, 0, 0), ir 2 = (n, it, 0, 0), . . . Thus, the staggered Dirac operator has non-diagonal terms 
that couple together states at different corners of the Brillouin zone. 
In momentum space, we can write the fermion matrix as 



p> 

M^ p » tP > = y^-2zsin(p^) 6 p » +7Tli ,p> +2mS p n 



(16) 
(17) 
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Then, we have for M 

MK".v'Mr 



4 ^ X! sin(p M ) sm(p") 8 p » + 4:m 2 S Ptpl/ 

4^^sin(p M )sin(p i , + (tt^)^ - (7r„)„) + 4ra 5. 

4^ X sin(p M ) sin(p v ) V',p+*>-f* + 
4^sin 2 (p At ) ( 5 p 



(18) 
(19) 



(20) 



Noticing that p + — ~k v = p + "K v — tt^ (mod 2tt), and interchanging /i and v labels in the second piece of Eq. (|20[) . 
we see that all the off-diagonal pieces of M} M. cancel, and we are left only with the diagonal piece, 



p" ■ 



(21) 



As expected, we have a hard lower bound on the eigenvalue spectrum of )? min = 4m 2 . The upper bound y? max = 
16 + 4m 2 is realized when p = (| , ~ , | , ~). 

The situation for free-field p4 fermions is similar, but a little bit more complicated. Here, the p4 Dirac operator is 
given by 

Mx{x) = 2mx(x) + ci, 53^(a:) 0d> + n) - \(x - m)) + ( 22 ) 
+Cl - 2 ^2 ^2 ''mW + M + + x(a; + // - 2z/) ~x{x- fi + 2v) - x(x - fi - 2u)) , 

where ci.o = § and c% t 2 — -kj- Again, we can examine how M. acts on the momentum eigenstate in Eq. (114|) . As 
before, we see that we have off-diagonal pieces that come about as a direct result of the presence of staggered phases. 



Or, in matrix form, 



M X (x) = J2x(p)e lip+ ^ } - x ih^p) + 2m X (p)e^ x , 

v- 

h^ip) = 2ci,q sin(p M ) + 4ci,a / ] sin(?y) cos(2p„). 



(23) 
(24) 

(25) 



Calculating Ai'M., we see that the off-diagonal pieces are eliminated in the same way as in the naive staggered case. 
Thus, we are left with only diagonal terms, 



{M^M) pl ^ p = h l(p) <W" + 4m 2 6 p y, 



(26) 



We see that ~>? min = 4m 2 . Finding X 2 lax requires us to maximize the function Doing this, we find the 

maximum when two of the components of p are equal to 7r/2 and the other two components are equal to 0. For 
example, p ma x — {n/2, 0, tt/2, 0). This yields A 2 ,^ = 50/9 + 4m 2 for the p4 action. A similar calculation for the Naik 
action yields A 2 71QX = 196/9 + 4m 2 . 

For completeness, we also quote the eigenvectors of the free p4 Dirac operator. Using a slightly different method 
we find, 



X (x) = Xp(X) 



'■^2 T pp'(p) s[n Pp ( 2c i,o + 4ci, 2 X cos2p l/ J + (\-2m)8 pp i 



u° p , e^ x , 



(27) 
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where we implicitly sum over p' . Here we have used hypercube coordinates X and p, where X labels the hypercube 
and p is the offset within the hypercube, x = 2X + p with p^ = 0, 1. Furthermore, T pp , (j>) is defined by T pp , (p) := 



free p4 Dirac operator is 



]e tp ( p p ' and vP is a constant vector depending only on p. Finally, the eigenvalue of the 



A = 2m ± ' 



\ 



^2 sin 2 p^ 



2ci, 



4ci,2 cos 2p„ 



(28) 
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